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[l] The combination of the Wang-Sheeley-Arge (WSA) coronal model, ENLIL heliospherical model 
version 2.7, and CONED Model version 1.3 (WSA-ENLIL with CONED Model) was employed to form 
ensemble forecasts for 15 halo coronal mass ejections (halo CMEs). The input parameter distributions 
were formed from 100 sets of CME cone parameters derived from the CONED Model. The CONED 
Model used image processing along with the bootstrap approach to automatically calculate cone 
parameter distributions from SOHO/LASCO imagery based on techniques described by Pulkkinen et al. 
(2010). The input parameter distributions were used as input to WSA-ENLIL to calculate the temporal 
evolution of the CMEs, which were analyzed to determine the propagation times to the Li Lagrangian 
point and the maximum K p indices due to the impact of the CMEs on the Earth's magnetosphere. The 
Newell et al. (2007) K p index formula was employed to calculate the maximum K p indices based on the 
predicted solar wind parameters near Earth assuming two magnetic field orientations: a completely 
southward magnetic field and a uniformly distributed clock-angle in the Newell et al. (2007) K p index 
formula. The forecasts for 5 of the 15 events had accuracy such that the actual propagation time was 
within the ensemble average plus or minus one standard deviation. Using the completely southward 
magnetic field assumption, 10 of the 15 events contained the actual maximum K p index within the range 
of the ensemble forecast, compared to 9 of the 15 events when using a uniformly distributed clock angle. 
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forecasting of coronal mass ejections using the WSA-ENLIL with CONED Model, Space Weather, 11, 95-106, 
doi:10.1002/swe.20019. 


1. Introduction 

[ 2 ] Coronal mass ejections (CMEs) are the cause of the 
most severe geomagnetic storms [Gosling, 1993]. Geomag- 
netic storms can cause a variety of complications at Earth 
including degradation of satellite performance [Afraimovich 
et al, 2003] and disruption of electrical systems on the 
Earth's surface [Boteler et al, 1998]. Therefore, the space 
weather community has a great interest in predicting the 
arrival times and impacts of CMEs at Earth. 

[ 3 ] A number of models have been developed to model 
CME propagation. Some of the earlier models include 
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shock propagation models based on type II meter wave 
burst measurements such as the Shock Time of Arrival 
(STOA) model [Dryer, 1974] and the Interplanetary Shock 
Propagation Model (ISPM) [Smith and Dryer, 1990], and 
kinematic models based on fast and slow solar wind 
stream interaction such as the Hakamada-Akasofu-Fry 
(HAF) model [Hakamada and Akasofu, 1982; Fry et al, 2001]. 
Recently, empirical forecast models have been developed, 
including the model developed by Gopalswamy et al [2001] 
which treats the CME as a kinematic object experiencing 
accelerations or decelerations to match the ambient solar 
wind speed at distances near 1 AU. 

[ 4 ] The most current and advanced method of forecast- 
ing CMEs is based on numerically approximating the 
solutions to the magnetohydrodynamic (MHD) equations 
governing the motion of the CME over time. ENLIL is a 
time-dependent three-dimensional model which solves 
the MHD equations for plasma mass, momentum, mag- 
netic field, and energy density using a finite-difference 
approximation [Odstrcil and Pizzo, 1999], ENLIL can accept 
the output of the Wang-Sheeley-Arge (WSA) coronal 
model to form the inner boundary condition in the finite- 
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difference computations, which calculates the background 
solar wind solution and interplanetary magnetic field 
(IMF) based on remote solar magnetogram measurements 
[Arge and Pizzo, 2000]. ENLIL can also accept the output of 
the Cone Model to initialize the CME speed, angular 
width, and axis of propagation. 

[ 5 ] The Cone Model, developed by Zhao et al. [2002], 
assumes that the CME has the extent of a cone with con- 
stant angular width, propagates in a radial direction, and 
experiences self-similar expansion. A practical technique 
to manually determine the cone parameters from SOHO/ 
LASCO imagery was developed by Xie et al. [2004], Previ- 
ous analyses have been completed using the analytic Cone 
Model along with WSA-ENLIL to forecast the impacts and 
propagation times of CMEs to Earth [Odstrcil et al., 2004], 
and have showcased the effectiveness of the WSA-ENLIL 
with Cone Model combination [e.g., Taktakishvili et al, 
2009; Taktakishvili et al., 2010]. 

[6] The analytic Cone Model relies on a manual deter- 
mination of the CME outer extent from LASCO imagery. 
The development of the CONED Model, an automated 
version of the Cone Model, removed the user from the 
subjective process of manually determining the CME 
outer extent [Pulkkinen et al, 2010]. The CONED Model 
uses image processing to automatically determine the lo- 
cation of the CME mass from a time-series of LASCO 
images and then calculates a distribution of possible cone 
parameters using the bootstrap approach [Pulkkinen et al, 
2010]. The distribution of cone parameters allows for a 
dynamic quantification of the uncertainty of the cone 
parameters based on LASCO imagery, which will vary for 
each event. 

[ 7 ] The performance of the WSA-ENLIL with CONED 
Model has been analyzed with the median values of the 
cone parameter distributions used as input for a single 
WSA-ENLIL run. The Taktakishvili et al. [2011] analysis 
showed that the analytic Cone Model and CONED Model 
version 1.0 (automatic Cone Model) had reasonable 
agreement in the forecasts with a mean absolute propa- 
gation time forecast error at Earth of 6.9 h for the analytic 
Cone Model and 11.2 h for the CONED Model. The results 
of this analysis indicated that CONED Model version 1.0 
required improvements to the cone parameter estimations 
to match the accuracy of the analytic Cone Model. The 
performance of the WSA-ENLIL with CONED Model 
version 1.2 was analyzed by Palkenberg et al [2011], with 
the conclusion that the CME speed and angular width 
were underestimated by CONED Model version 1.2. 
CONED Model version 1.3 is the most current version of 
the CONED Model and has included a modification in the 
optimization routine to increase the CME speed and width 
estimations following the results of the Palkenberg et al. 
[2011] analysis. 

[8] With the production of the cone parameter distribu- 
tions from the CONED Model readily available, an en- 
semble forecast can be applied. The weather community 
has long known of the improvement in forecast accuracy 
due to the use of ensemble forecasting [Leith, 1974]. 


Ensemble forecasting also allows for a quantification of 
forecast uncertainty based on uncertainty in the mea- 
surements of the initial conditions, which is impossible for 
single forecasts. This quantification of forecast uncertainty 
could provide useful information to operational forecasts 
of CMEs. 

[ 9 ] We applied the ensemble forecasting technique to 
15 halo CMEs using the WSA-ENLIL with CONED Model 
version 1.3. The ensembles were created from 100 sets of 
initial states (cone parameters) derived from the CONED 
Model, which were used as input to WSA-ENLIL version 2.7 
to obtain distributions of propagation times to the L, 
Lagrangian point and distributions of maximum K p indices 
due to the impact of the CME on the Earth's magnetosphere. 

2. Brief Description of WSA, ENLIL, 
and CONED Model 

[ 10 ] In 2009, Pulkkinen et al. created the CONED Model, a 
tool which determines the cone parameters of CMEs from 
time series of LASCO C3 difference images automatically. 
The cone parameters are composed of the radial speed of 
the cone-front, the angular half width of the cone, the 
propagation axis (direction of propagation) of the cone, and 
the position of the CME leading edge as a function of time. 
The CONED Model uses image processing to automati- 
cally determine the shape of the CME mass from LASCO 
C3 difference images by filtering the images based on a 
brightness threshold, which differentiates the brighter 
CME mass location from the darker image background. 
The bootstrap approach is used to determine the confi- 
dence intervals for the calculated cone parameters by 
randomly selecting 300 points from the CME mass loca- 
tions in the filtered LASCO images and then calculating 
the best-fitting cone to the randomly selected points. The 
process can be repeated any number of times, to create a 
distribution of cone parameters for the particular CME. 

[ 11 ] The Wang-Sheely-Arge (WSA) model is an empiri- 
cal model used to calculate background solar wind speed 
and IMF [Arge and Pizzo, 2000]. The model calculates the 
magnetic field between the solar surface and a boundary 
sphere where the magnetic field is assumed to be radial, 
based on synoptic magnetogram data. The solar wind 
speed at the boundary sphere is calculated using an 
inverse relationship between the solar wind speed and the 
magnetic expansion factor. WSA provides synoptic maps 
of the radial magnetic field and speed that are used by 
ENLIL to calculate MHD boundary conditions at 21.5 solar 
radii. 

[ 12 ] After the input parameters are obtained from the 
CONED Model and the synoptic maps are obtained from 
the WSA model, ENLIL approximates the time-dependent 
solution to the MHD equations governing the plasma from 
21.5 solar radii to an appropriate outer boundary (1.1 AU 
for quickly analyzing the effects of a CME near Earth). 
ENLIL utilizes a modified Total-Variation-Diminishing 
Lax Friedrichs (TVDLF) finite difference scheme to ap- 
proximate the solution to the partial differential MHD 
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equations [Toth and OAstrcil, 1996]. The current version of 
ENLIL assumes no internal magnetic field structure to the 
CME but allows for the distortion of the IMF by shock 
deflection (if a shock develops), draping around ejecta and 
stretching in the rarefaction region. 

3. Ensemble Forecasting 

[ 13 ] Formally, ensemble forecasting can be described by 
a transition from a probability distribution of initial states, 
p(vt\dt), given a set of observations, o t , to a probability 
distribution of future states, p(v t+z \b t ): 

p(v t +z\o t ) = j r(v t+z \v t )p(v t \o t )dv t , (1) 

where v, is the initial state, v t+z is the future state, and 
r(v f+ r |tif) is the transition probability associated with the 
forecasting model [DelSole, 2005]. For a deterministic 
model (a model which provides the same result if run mul- 
tiple times with the same set of initial conditions), such as 
ENLIL, the transition probability can be described by a 
delta function. 

[ 14 ] The ensemble forecast distribution for a particular 
set of observations can be calculated from the probability 
distribution of future states. The ensemble forecast distri- 
bution can be statistically analyzed to obtain the most 
likely value of a particular parameter, along with the as- 
sociated uncertainty of the value. The range of the en- 
semble forecast distribution provides the range of possible 
outcomes for a given set of observations. 

[ 15 ] For this analysis, an ensemble forecast was calcu- 
lated for a total of 15 halo CMEs using the WSA-ENLIL 
version 2.7 with CONED Model version 1.3. For each 
CME, the CONED Model was used to sample 100 sets of 
initial conditions from the probability distribution of initial 
states based on three LASCO C3 difference images of the 
CME eruption with a temporal spread of at least 1 h be- 
tween the first and last images. The 100 sets of initial 
conditions were then used as input to WSA-ENLIL to 
obtain the probability distributions of future states, which 
were used to calculate the ensemble forecast distributions. 
The ensemble forecasting process used for this analysis 
can be summarized by 

LASCO C3 Difference Images —> dt — > Coned Model 
-> p(vt\o t ) -> WSA - ENLIL -» p(v t+z \b t ) 

— > Ensemble Forecast Distribution. (2) 

[ 16 ] Two parameters were calculated from p(v t+r \o t ) to 
obtain the ensemble forecast distributions: the propaga- 
tion time of the CME to the Li Lagrangian point and the 
maximum K p index due to the CME impact on the Earth's 
magnetosphere. For this analysis, the resolution of the 
computational grid used by ENLIL (r x 0 x cj> = 160 x 30 x 90) 
placed the L, Lagrangian point and Earth in the same grid 
cell, so the computed propagation time to Earth was the 
same as the computed propagation time to the L a Lagrang- 
ian point. 


[ 17 ] The arrival time of the CME at Earth was selected to 
be the time at which the first temporal derivative of the 
solar wind dynamic pressure, calculated from ENLIL out- 
put, experienced a large increase in magnitude. The ar- 
rival time could also be considered to be the time at which 
the second temporal derivative of the dynamic pressure 
was a maximum. The arrival times computed using the 
first and second derivatives were compared to ensure no 
falsely triggered arrival times. Although not shown, the 
results were in good agreement. 

[is] The maximum K p indices were found using the 
Newell et al. [2007] K p formula: 

K p = 0.0002947 u 4 / 3 B P /3 sin 8 / 3 (e c /2) + 1, (3) 

where v is the speed (km/s) at which the IMF lines ap- 
proach the magnetopause and can be approximated by 
the solar wind speed, B T is the magnitude of the IMF (nT), 
and 0, is the IMF clock angle defined by 0 C = arctan(B y IB z ). 
The magnetic field orientation of the CME "cloud" is not 
computed by ENLIL, so two magnetic field orientations 
were assumed: First, the magnetic field was assumed to 
be completely southward ( 9 C = n), in order to calculate the 
worst-case scenario. Second, the expected value of the 
clock-angle term in the Newell et al. [2007] formula was 
used, assuming the clock-angle was randomly oriented 
with a uniform distribution. For a uniformly distributed 
clock-angle, the expected value of the clock-angle term 
was calculated to be 



[ 19 ] As a first-order estimation, the instantaneous (point 
by point) K p was calculated from the WSA-ENLIL with 
CONED Model data instead of using the 6h integration 
time described by Newell et al. [2007]. With no magnetic 
field orientation information available from ENLIL, the 
instantaneous K p values were found to be nearly identical 
to the K p values calculated using the 6 h integration time 
(results not shown). The K p indices were rounded to the 
nearest integer value, and any calculations of the maxi- 
mum K p index using the Newell et al. [2007] formula ex- 
ceeding nine were truncated to nine. 

[ 20 ] Eight of the CMEs studied in this analysis were se- 
lected from the Taktakishvili et al. [2011] analysis, based on 
CMEs which caused particularly large geomagnetic 
storms. The other seven CMEs had associated maximum 
K p values less than eight and had no other halo CMEs 
within plus or minus 2 days from the eruption day of the 
CME. All of the CMEs were required to have clear LASCO 
C3 difference images to run the CONED Model and clear 
ACE data to determine the actual arrival times of the 
CMEs at the L, Lagrangian point. The CMEs were also 
selected with a large variety of associated solar flare loca- 
tions in order to analyze the performance of the model 
with CMEs initiated from a variety of solar locations. 
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Table 1. The Start Date and Times, Actual Propagation Times as Measured by ACE, Actual Maximum K p Indices, and the Loca- 
tions of the Associated Solar Flares for the 15 CMEs Analyzed. The CMEs Were Also Labeled With an Event Number for Easy 
Reference 


Event 

Number 

CME Start Date 
(Date/Month/Year) 

CME Start 
Time (UT) 

Propagation Time 
to ACE (HH:MM) 

Maximum K p 

Associated Solar 
Flare Location 

1 

03 May 1999 

06:06 

56:50 

3 

N15E32 

2 

04 April 2000 

16:32 

47:30 

9 

N16W66 

3 

14 July 2000 

10:54 

27:20 

9 

N22W07 

4 

29 March 2001 

10:26 

37:50 

9 

N20W19 

5 

10 April 2001 

05:30 

33:50 

8 

S23W09 

6 

24 September 2001 

10:30 

33:30 

7 

S16E23 

7 

09 October 2001 

11:30 

52:45 

6 

S28E08 

8 

04 November 2001 

16:35 

32:40 

9 

N06W18 

9 

17 November 2001 

05:30 

60:00 

4 

S13E42 

10 

28 October 2003 

11:30 

18:20 

9 

S16E08 

11 

29 October 2003 

20:54 

19:50 

9 

S15W02 

12 

20 July 2004 

13:31 

44:20 

7 

N10E35 

13 

06 November 2004 

02:06 

39:40 

9 

N07E00 

14 

03 December 2004 

00:26 

54:20 

4 

N09E03 

15 

03 April 2010 

10:34 

45:15 

8 

S25E00 


[ 21 ] The calculated propagation times were compared to 
the actual propagation times derived from ACE measure- 
ments, with 10 min resolution. The actual arrival times 
calculated directly from ACE data were compared to the 
arrival times logged in the National Oceanic and Atmo- 
spheric Administration (NOAA) Space Weather Predic- 
tion Center's (SWPC's) historical weekly reports (http:// 
www.swpc.noaa.gov/ftpmenu/warehouse.html) to ensure 
consistency. 

[ 22 ] The calculated maximum K p indices were compared 
to the actual ground-based maximum K p values, with inte- 
ger resolution. The actual maximum K p indices were found 
using NASA's OMNIWeb database (http://omniweb.gsfc. 
nasa.gov/form/dxl.html). The solar flare locations, derived 
from USAF/SOON observations, were used to approximate 
the location for the CME eruptions. The measured values 
for the actual propagation times, maximum K p indices, and 
locations of the associated solar flares are displayed in 
Table 1. 

4. Results 

4.1. Single Event 

[ 23 ] The ensemble distributions for event 4 (29 March 
2001 CME) are displayed in this section as an example of 
the results obtained for each of the CMEs analyzed. The 
LASCO/C3 images used to create the cone parameter 
distributions are displayed in Figure 1, and the cone pa- 
rameter distributions are displayed in Figure 2 and 
Table 2. 

[ 24 ] Each of the 100 sets of cone parameters mapped to a 
separate forecast (Figure 3). The set of 100 forecasts created 
by the 100 sets of cone parameters formed the ensemble 
forecast (Figure 4). The ensemble forecast was analyzed 
to obtain a single forecast along with the uncertainty in 
the forecast, which was calculated from the width of the 
ensemble (Table 2). The propagation time ensemble aver- 
age was used as the single forecast for the propagation 


time, and the rounded average (to the nearest integer) of 
the maximum K p index ensemble was used as the single 
forecast for the maximum K p index. The standard devia- 
tions and ranges were used to describe the width of the 
ensembles. 

[ 25 ] For this event, the average propagation time was 
36.4 h with a standard deviation of 5.8 h. With an actual 
propagation time of 37.8 h, the absolute forecast error 
(absolute value of the ensemble average minus the actual 
value) was 1.4 h. Using the completely southward mag- 
netic field assumption, the average maximum K p was cal- 
culated to be 9, with a range of 2 (the standard deviation 
was 0.2, which was too small to represent the uncertainty 
for the integer resolution used for the K p ). Using a uni- 
formly distributed clock angle, the ensemble average was 
7, with a range of 4. With an actual maximum K p index of 
9, the completely southward magnetic field predicted the 
magnitude perfectly, while using a uniformly distributed 
clock angle underestimated the magnitude by 2. 

4.2. Cone Parameters 

[ 26 ] The cone parameter distributions for each of the 15 
CMEs are displayed in Figure 5. The 15 cone parameter 
distributions were averaged to obtain an overall measure 
of the uncertainty derived from LASCO imagery via the 
CONED Model (Table 3). 

[ 27 ] The averages of the CONED Model speed distribu- 
tions displayed a strong positive correlation with the 
LASCO first-order plane of sky (POS) speeds (Figure 6). 
The LASCO first-order POS speeds were the first-order 
(linear) fits to the leading edge of the CMEs in the LASCO 
POS imagery (http://cdaw.gsfc.nasa.gov/CME_list/) and 
provided an estimate of the two-dimensional POS speed 
of the CMEs. The strong correlation was expected because 
of the connection between the projected (POS) speed and 
the radial speed, and the fact that this analysis focused on 
geo-effective CMEs with a limited range of latitudes and 
longitudes. The average ratio of the CONED Model 
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Figure 1. The LASCO C3 images used to calculate the cone parameter distributions for 
event 4 (29 March 2001). The format of the image time stamps follows: YYYYMMDD_ 
HHMM. Note the 1 h time elapse from the first and last images used to calculate the cone 
parameter distributions. 


average speeds to the LASCO first-order POS speeds was 
1.2, which was consistent with the fact that the CONED 
Model predicted the radial (three-dimensional) speed, 
while the LASCO first-order POS speed was a measure of 
the two-dimensional POS speed. Additionally, the average 
ratio of the CONED Model radial speeds to the CONED 
Model POS speeds was calculated to be 1.1, which agreed 
well with the ratio of the CONED Model radial speeds to 
the LASCO first-order POS speeds. 

[ 28 ] While the location of the associated solar flare is not 
necessarily an indicator of the direction of the CME 
propagation, the CONED Model tended to push the 
propagation axes of the CMEs towards the Sun-Earth 
line relative to the locations of the associated solar flares 
(Figure 7). The CONED Model calculated propagation 
axes with average latitude or longitude magnitudes 
greater than 10° for 4 of the 15 CMEs, while 13 of the 15 
associated solar flare location latitude or longitude mag- 
nitudes were greater than 10°. CONED Model version 1.3 


does not use the solar flare location as a constraint in the 
cone parameter estimations, but the results of this analysis 
have shown the necessity to use additional information in 
the calculations (efforts are currently underway to add this 
capability). 

4.3. Propagation Time 

[ 29 ] The ensemble forecasts contained 5 of the 15 actual 
propagation times within one standard deviation of the av- 
erage (Figure 8). These five CMEs had actual propagation 
times between 30 and 46 h. There were two other CMEs with 
similar propagation times that were not forecast within one 
standard deviation. The actual propagation time for 8 of the 
15 ensemble forecasts were within of the range of the en- 
semble distribution (Figure 9). Of the eight forecasts, seven 
were for CMEs with actual propagation times between 30 
and 46 h, and the remaining forecast was for a CME with an 
actual propagation time around 53 h. 
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Figure 2. The cone parameter distributions, derived from CONED Model version 1.3, for 
event 4 (29 March 2001 CME). The count number is the number of times a cone parameter 
was within the range of a bin, for the 100 sets of cone parameters. 


[ 30 ] The average of the ensemble standard deviations 
was calculated to be 4.6 ± 0.7 h (Table 3). This quantifica- 
tion of the propagation time uncertainty was based strictly 
on the uncertainty in the initial conditions calculated from 
LASCO imagery of the particular CMEs. Another measure 
of the propagation time uncertainty was the range of the 
ensembles, which were averaged to obtain 22.2 ± 5.2 h. 
While this was too large of an uncertainty to be useful for 
operational forecasts, it was an important metric to ana- 
lyze the overall performance of the ensemble forecasting 
technique. The quantification of the uncertainty provided 
an estimate of the widths of the ensembles and was not a 
measure of the forecast error. 

[ 31 ] The ensemble forecasts for the eight events with 
actual propagation times between 27 h and 46 h were the 
most accurate of the 15 forecasts, with all of the absolute 
forecast errors less than 8 h (Figure 10). The forecast errors 
for the two fast CMEs, events 10 and 11 (28 October 2003 
and 29 October 2003 CMEs), were around 9h indicating 
that the CME speeds were underestimated. The absolute 
forecast errors for events with actual propagation times 
less than 46 h were all less than 10 h. 

[ 32 ] The forecast errors for the slow CMEs, with actual 
propagation times over 46 h, were all less than — 10 h. This 
indicated that the ensemble forecasts greatly under- 
estimated the propagation times of the slower CMEs. For 
the events with actual propagation times greater than 50 h, 
the absolute forecast error increased as the actual propaga- 
tion time increased. The slowest event (event 9) had an 
actual propagation time of 60 h and a forecast error of -25 h. 


The large forecasting errors for the slower CMEs were due 
the fact that CONED Model 1.3 was tuned to more accu- 
rately forecast faster (more significant from the space 
weather point of view) events following the Falkenberg et al. 
[2011] study, forcing the CONED Model to select faster 
speeds. 

[ 33 ] The mean absolute forecast error, for the 15 CMEs, 
was calculated to be 9.1 ± 7.1 h (Table 3). While the set of 
events between this analysis and the Taktakishvili et al. 
[2011] analysis were not identical, the ensemble forecast- 
ing mean absolute forecast error was greater than the 
mean absolute error of 6.9 ± 4.2 h found by Taktakishvili 
et al. [2011] using single ENLIL runs with the analytical 
Cone Model but was less than the 11.2 ± 7.2 h mean abso- 
lute error found by Taktakishvili et al. [2011] using single 
ENLIL runs with the median values of the cone parameters 
derived from CONED Model version 1.0. It must be noted 
that the Taktakishvili et al. [2011] analysis focused on CMEs 
associated with particularly large geomagnetic storms, 
while this analysis included multiple CMEs that produced 
relatively weak geomagnetic storms (K p < 4), which could 
be considered less significant from the space weather point 
of view. Therefore, the forecast errors may not be directly 
comparable between the analyses. 

4.4. Maximum K p Index 

[ 34 ] The ensemble forecast tended to overestimate the 
magnitude of the impacts of the CMEs by forecasting a 
maximum K p of 9 for all 15 CMEs using the completely 
southward magnetic field assumption (Figure 11). Of the 
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Table 2. The Average and Standard Deviation Time, and Maximum K p Distributions, for Event 4 (29 March 2001 CME) 


Parameter 

Average 

Standard Deviation 

Minimum 

Maximum 

Range 

Speed (km/s) 

1444.3 

304.9 

848.0 

2256.0 

1408.0 

Angular Half Width (deg) 

53.8 

10.3 

33.0 

82.0 

49.0 

Latitude (deg) 

-0.1 

0.3 

-1.0 

0.0 

1.0 

Longitude (deg) 

0.0 

0.2 

0.0 

1.0 

1.0 

Propagation time (hours) 

36.4 

5.8 

26.8 

57.3 

30.5 

Southward IMF K p 

9 

- 

7 

9 

2 

Uniform clock angle K p 

7 

- 

4 

8 

4 


15 CMEs, seven had an actual maximum K p of 9, and three 
had an actual maximum K p less than 5. The ensemble 
forecasts for 10 of the 15 CMEs had accuracy such that the 
actual maximum K p was within the range of the ensemble. 
Of the eight events with actual maximum K p indices less 
than 9, three forecasts contained the actual maximum K p 
inside of the range of the ensemble. 

[ 35 ] The average of the ensemble ranges, for the 15 
events, was calculated to be 0.7 ±1.3 (Table 3). The en- 
semble ranges were zero for all but four events due to the 
overestimation of the maximum K p values and the fact that 
any maximum K p calculation exceeding 9 was truncated to 
9. This provided a quantification of the uncertainty in the 
maximum K p calculations and was not a measure of the 
forecast error. The mean absolute forecast error, for all 
15 events, was calculated to be 1.7 ± 2.1. 

[ 36 ] The completely southward magnetic field maximum 
K p forecast is displayed along with the propagation time, 
per event, in Figure 12. The events with the largest prop- 
agation time errors also had the largest maximum K p 
errors. This was due to overestimations of the CME speeds 


for these particular events, which forecast the arrival times 
too early, and the maximum K p indices too large. 

[ 37 ] The maximum K p forecasts for 8 of the 15 events 
were lowered by using a uniformly distributed clock angle, 
relative to the completely southward magnetic field 
forecasts (Figure 11). The forecasts for 9 of the 15 events had 
accuracy such that the actual maximum K p was within the 
range of the ensemble. The forecasts for 4 of the 8 events 
with actual maximum K p indices less than 9 had the actual 
maximum K p within the range of the ensemble, which was 
slightly better than the 3 of 8 for the completely southward 
magnetic field forecasts. But, the forecasts using a uni- 
formly distributed clock-angle tended to underestimate the 
maximum K p indices for the events with actual maximum 
K p indices of 9. The average of the ensemble ranges, for the 
15 CMEs, was calculated to be 1.5 ± 1.4 (Table 3). The range 
was 0 for 4 of the events, compared to 11 events for the 
completely southward magnetic field forecasts. 

[38] The mean absolute forecast error for the 15 CMEs 
using a uniformly distributed clock angle was calculated to 
be 1.8, which was slightly larger than the 1.7 calculated 
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Figure 3. The cone parameters and propagation time forecasts for each of the 100 sets of 
parameters composing the ensemble for event 4 (29 March 2001 CME). Each set of cone para- 
meters maps to a propagation time, which shows the relationship between the cone para- 
meters and the ensemble forecast. 
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Figure 4. The maximum K p index and propagation 
time ensemble forecasts, for event 4 (29 March 2001 
CME). The count number is the number of times a cone 
parameter was within the range of a bin, for the 100 sets 
of forecasts derived from the 100 sets of cone para- 
meters. The "southward IMF" describes the maximum 
K p index calculations using a completely southward 
magnetic field, and the "uniformly distributed clock an- 
gle" describes the calculations using a uniformly dis- 
tributed clock angle in the Newell et al. [2007] formula. 

using the completely southward magnetic field forecasts. 
The mean absolute forecast error for the events with actual 
maximum K p indices of 9 was 0.7, which was greater than 


the 0.0 for the completely southward magnetic field fore- 
casts. The mean absolute forecast error for the events with 
actual maximum K p indices less than 9 was 2.8, which was 
less than the 3.1 for the completely southward magnetic 
field forecasts. This indicated that the forecasts completed 
using a uniformly distributed clock angle were less accu- 
rate than the forecasts completed using a completely 
southward magnetic field for the events with actual max- 
imum K p indices of 9 but were more accurate for the 
events with actual maximum K p indices less than 9. This 
follows from the fact that the magnetic field orientation of 
the CMEs that produced the weak geomagnetic storms 
were not conducive to producing large geomagnetic 
storms (there was not a large southward component), and 
the assumption that the magnetic field was completely 
southward tended to overestimate the magnitude of the K p 
for these events. 

[ 39 ] Using ACE data in the Newell et al. [2007] formula, 
along with the 6h integration time described by Newell 
et al. [2007] (with SOHO Proton Monitor data to estimate 
the solar wind speed for events 8, 10, and 11, which had 
missing solar wind speed information from ACE due to 
proton events), the mean absolute error for the 15 events 
was calculated to be 0.5. This indicated that the Newell et al. 
[2007] formula used with actual solar wind data (including 
magnetic field orientation) as input could produce accurate 
K p estimates. Therefore, the main source of the errors in the 
maximum K p estimations from WSA-ENLIL with CONED 
Model were due to errors in the estimations of the solar 






Figure 5. The average and standard deviation of the cone parameter distributions, derived 
from CONED Model version 1.3, for each event. The CONED Model used three LASCO C3 
images of the CME eruption to develop the distributions of cone parameters, for each CME. 


102 


EMMONS ET AL.: CME ENSEMBLE FORECASTING 


Table 3. The Averages of the Standard Deviations, Ranges, and Forecast Errors for the Cone Parameter, Propagation Time, and 
Maximum K p Distributions, Over All 15 Event 


Parameter 

Average of the Standard Deviations 

Average of the Ranges 

Average of the Forecast Errors 

Speed (km/s) 

323.6 

1568.3 

- 

Angular Half Width (deg) 

9.1 

43.3 

- 

Latitude (deg) 

1.0 

5.3 

- 

Longitude (deg) 

1.3 

6.8 

- 

Propagation time (hours) 

4.6 

22.2 

9.1 

Southward IMF K p 

- 

0.7 

1.7 

Uniform clock angle K p 

- 

1.5 

1.8 


wind parameters near Earth and the fact that no magnetic 
field orientation information was available. 


5. Conclusions 

[ 40 ] Ensemble simulations of CME impact and propa- 
gation to Earth were produced for 15 halo CMEs using the 
WSA-ENLIL version 2.7 with CONED Model version 1.3, 
helping to close the gap between terrestrial and space 
weather forecasting. The ensemble forecasts consisted of 
the propagation times to the Li Lagrangian point and the 
associated maximum K p indices due to the impact of the 
CMEs on the Earth's magnetosphere. The ensemble fore- 
casts were formed by using 100 sets of cone parameters, 
derived from LASCO C3 imagery via the CONED Model, 
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Figure 6. The CONED Model average speeds versus 
the LASCO first-order POS speeds, with the event 
number as the label and the standard deviations of 
the ensembles as the error bars. The two methods of 
calculating the CME speeds were strongly correlated, 
with the average ratio of the CONED Model average 
speeds to the LASCO first-order POS speeds calculated 
to be 1.2. 


corr=0.90, p-val=0.00 
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as input to WSA-ENLIL to produce propagation time and 
maximum K p index distributions. 

[ 41 ] The propagation time ensemble forecasts estimated 
5 of 15 events with accuracy such that the actual propa- 
gation time was within the ensemble average plus or mi- 
nus the ensemble standard deviation. The forecasts for 8 of 
15 events had accuracy such that the actual propagation 
time was within the range of the ensemble. The mean 
absolute forecast error, for the 15 CMEs, was calculated to 
be 9.1 h. The majority of the error came from the slow 
events because CONED Model 1.3 was tuned to more ac- 
curately forecast faster events following the Falkenberg 
et al. [2011] study. 

[ 42 ] Perhaps the most important result of this analysis 
was the dynamic quantification of the forecast uncertainty 
derived strictly from measurements (LASCO imagery) of 
particular CMEs. The average of the propagation time 
standard deviations, over all 15 events, was calculated to 
be 4.6 h, and the average of the ensemble ranges was 


Corr=0.70, p-val=0.00 



Figure 7. The CONED Model average longitudes and 
latitudes along with the associated solar flare latitudes 
and longitudes, with the event numbers as the labels 
and the standard deviations of the ensembles as the 
error bars. Relative to the associated solar flare loca- 
tions, the CONED Model propagation axes tended to 
be pushed towards the Sun-Earth line. 
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Figure 8. The averages and standard deviations of the 
propagation time ensembles versus the actual propaga- 
tion times, with the event number as the label. The 
actual propagation time was within the average plus or 
minus one standard deviation for 5 of the 15 events. 


Figure 10. The propagation time forecast error versus 
the actual propagation time, with the error bars as one 
standard deviation and the labels as the event number. 
The forecasts for the events with actual propagation 
times between 27 and 46 h were the most accurate, 
while the largest forecast errors were for the events 
with actual propagation times greater than 46 h. 


calculated to be 22.2 h. While these values were not a 
measure of the forecast accuracy, they did provide a 
measure of the uncertainty in the forecasts based on the 
uncertainty in the measurements of the initial conditions, 
which would be useful information in operational fore- 
casts of CMEs. 

[ 43 ] The maximum K p indices were calculated using the 
maximum K p index formula created by Newell et al. [2007] 
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Figure 9. The averages and ranges of the propagation 
time ensembles versus the actual propagation times, 
with the event number as the label. The actual propaga- 
tion time was within the ensemble range for 8 of the 
15 events. 


using two different magnetic field orientations. Using the 
completely southward magnetic field, the ensemble fore- 
cast predicted maximum K p indices of 9 for all events, 
which was an overestimation for many of the events. The 
forecasts for 10 of the 15 events had accuracy such that the 
actual maximum K p index was within the range of the 
ensemble forecast. The mean absolute forecast error was 
calculated to be 1.7, and the average of the ensemble 
ranges was 0.7. 

[44] Using a uniformly distributed clock angle, 9 of the 
15 events were forecast with accuracy such that the actual 
maximum K p index was within the range of the ensemble. 
The mean absolute forecast error was calculated to be 1.8, 
and the average of the ensemble ranges was 1.5. The 
forecasts using a uniformly distributed clock angle were 
more accurate than the forecasts assuming the magnetic 
field was completely southward for the events with actual 
maximum K p indices less than 9 but were less accurate for 
the events with actual maximum K p indices of 9. 

[45] The next step in ensemble forecasting of CMEs is to 
advance the CONED Model to produce more accurate 
estimations of the initial CME parameters. Recently, Jacobs 
et al. [2012] developed more advanced segmentation tech- 
niques to use in the model. The current model develop- 
ment is focused on allowing the model to process multiple 
transients in the field of view and allowing the location of 
the CME eruption (associated solar flare location) as well as 
any propagation information derived from STEREO to be 
taken into account when calculating the cone parameters. 
Improvements in the estimations of the cone parameters 
should improve the overall accuracy of the ensemble 
forecasting technique. 
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Figure 11. The average and range of the maximum K p ensembles, per event, using a uni- 
formly distributed clock angle in the Newell et al. [2007] formula (top) and the assumption that 
the magnetic field is completely southward (bottom). The red squares with error bars are the 
ensemble forecasts, the bars are the actual maximum K p indices, and the blue triangles are 
the calculated maximum K p indices using ACE solar wind data in the Newell et al. [2007] 
formula. 


[46] Additionally, the K p forecasts will be recalculated 
using the actual magnetic field orientations of the CMEs 
along with the WSA-ENLIL with CONED Model data to 
determine the improvement in the forecast accuracy if the 


actual magnetic field information was available. This will 
help determine the usefulness of accurate clock-angle esti- 
mates for K p forecasts using the WSA-ENLIL with CONED 
Model. 




Event 

Figure 12. The average and range of the maximum K p index ensemble using the completely 
southward magnetic field assumption, along with the average and standard deviation of the 
propagation time, per event, for the 15 CMEs. The ensemble forecasts and uncertainties are 
the points with the error bars, and the bars are the actual values. The events with the largest 
propagation time errors also had the largest maximum K p errors, which was due to overesti- 
mations of the CME speeds. 
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